Semivariogram Subroutine

private subroutine Semivariogram(stations)

generate a Semivariogram from a given sample. The sample point pairs are ordered into even-width bin, separated by the euclidean distance of the point pairs. The Semivariance in the bin is calculated by the Matheron estimator. If number of lags and lag max distance are not given they are automatically computed or set to default value.

References:

Matheron, G. (1965) Les variables régionalisées et leur estimation: Une application de la théorie des fonctions aléatoires aux sciences de la nature. Masson, Paris.

Arguments

Type IntentOptional Attributes Name
type(ObservationalNetwork), intent(in) :: stations

Variables

Type Visibility Attributes Name Initial
real(kind=float), public, ALLOCATABLE :: distCopy(:,:)
real(kind=float), public, ALLOCATABLE :: empVariogramCopy(:,:)
integer(kind=short), public :: i
integer(kind=short), public :: index
integer(kind=short), public :: j
integer(kind=short), public :: k
integer(kind=short), public :: m
integer(kind=short), public, ALLOCATABLE :: pairNumberCopy(:,:)
real(kind=float), public :: zj

observation at point j, and k
real(kind=float), public :: zk

observation at point j, and k

Source Code

SUBROUTINE Semivariogram &
!
(stations)

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: stations

!Local declarations
INTEGER (KIND = short) :: i, j, k, m, index
REAL (KIND = float) :: zj, zk !!observation at point j, and k
REAL (KIND = float), ALLOCATABLE :: empVariogramCopy (:,:), distCopy (:,:) 
INTEGER (KIND = short), ALLOCATABLE :: pairNumberCopy (:,:)

!---------------------------end of declarations--------------------------------

!generate empirical semivariogram

DO j = 1, stations % countObs !loop through pairs
   DO k = 1, stations % countObs
      IF ( j == k ) CYCLE !skip same point pairs, distance is zero  
      IF ( k <  j ) CYCLE !this pair is already in the vector; pair(4,2) is the same as pair(2,4)
     
      DO m = 1, ndir
         DO i = 1, lagNumber (m) !loop through lags
                          
                IF ( ndir == 1 ) THEN !isotrophic analysis, direction is not considered
            
                    IF (dist_pairs (j,k) >= dist (m,i) - lagTolerance (m) .AND. &
                        dist_pairs (j,k) < dist (m,i) + lagTolerance (m) ) THEN
                 
                         zj = stations % obs (j) % value
                         zk = stations % obs (k) % value
                 
                         pairNumber (m,i) = pairNumber (m,i) + 1
                         empVariogram (m,i) = empVariogram (m,i) + ( zj - zk )**2.
                
                    END IF
                        
                ELSE !anisotrophy analysis, create bins considering direction
                    IF ( dist_pairs (j,k) >= dist (m,i) - lagTolerance (m) .AND. &
                            dist_pairs (j,k) < dist (m,i) + lagTolerance (m) ) THEN
                        
                        
                       IF (  ( m == 1 .AND. dir_pairs (j,k) < pi / 8. )   .OR. &
                          ( m == 1 .AND. dir_pairs (j,k) >= 15./8. * pi ) .OR. &
                          ( m == 2 .AND. dir_pairs (j,k) >= pi / 8. .AND. dir_pairs (j,k) < 3./8. * pi )    .OR. &
                          ( m == 3 .AND. dir_pairs (j,k) >= 3./8. * pi .AND. dir_pairs (j,k) <= 1./2. * pi ) .OR. &
                          ( m == 3 .AND. dir_pairs (j,k) > 3./2. * pi .AND. dir_pairs (j,k) < 13./8. * pi ) .OR. &
                          ( m == 4 .AND. dir_pairs (j,k) >= 13./8. * pi .AND. dir_pairs (j,k) < 15./8. * pi ) ) THEN
                        
                             zj = stations % obs (j) % value
                             zk = stations % obs (k) % value
                             pairNumber (m,i) = pairNumber (m,i) + 1
                             empVariogram (m,i) = empVariogram (m,i) + ( zj - zk )**2.
                    
                        END IF
                     END IF  
                        
                END IF 
            
         END DO
      ENDDO     
   END DO
END DO


DO m = 1, ndir !loop through directions
  DO i = 1, lagNumber (m) !loop through lags
     empVariogram (m,i) = empVariogram (m,i) / ( 2. * pairNumber (m,i) )
  END DO
END DO

RETURN

END SUBROUTINE Semivariogram